
clear
cd('Z:\Conformity Effects\Matlab')

% Set parameters
global p_values gammapar

% Load data into memory
data = csvread('zip3 data.csv');
data=data([1:184,186:end],:);
zip3=data(:,1); Ndem=data(:,2); Nrep=data(:,3); N=data(:,4); 
clear data;

% Maximum likelihood estimation
gammapar=0.11;
options = optimset('Display','iter','TolFun',1e-3,'TolX',1e-3);
initial=[6.0143 0.8754];
coef = fminsearch(@(x) lik(x,Ndem,Nrep,N),initial,options)

% Counterfactual simulations with estimated parameters
gammapar=0.11;
coef=[5.2 0.87]; % Estimated parameters
lik(coef,Ndem,Nrep,N) % Computes the auxiliary p_values
rng(123);
theta=coef(1);
alpha=coef(2);
sim_Mdem=round(betarnd(theta,theta,size(N)).*N);
sim_Mrep=N-sim_Mdem;
sim_Ndem=zeros(size(N)); sim_Nrep=sim_Ndem;
sim_pdem=p_values.*(sim_Mdem./sim_Mrep).^gammapar;
sim_prep=alpha.*p_values.*(sim_Mrep./sim_Mdem).^gammapar;
sim1_pdem=p_values;
sim1_prep=alpha.*p_values;
Niter=1000; aux_sim_Ndem=zeros(size(N,1),Niter); aux_sim_Nrep=aux_sim_Ndem; aux_sim1_Ndem=aux_sim_Ndem; aux_sim1_Nrep=aux_sim_Ndem;
for i=1:1:Niter
    aux_sim_Ndem(:,i)=randbinom(sim_Mdem',sim_pdem')';
    aux_sim_Nrep(:,i)=randbinom(sim_Mrep',sim_prep')';
    aux_sim1_Ndem(:,i)=randbinom(sim_Mdem',sim1_pdem')';
    aux_sim1_Nrep(:,i)=randbinom(sim_Mrep',sim1_prep')';
end
sim_Ndem=mean(aux_sim_Ndem,2); sim_Nrep=mean(aux_sim_Nrep,2); sim1_Ndem=mean(aux_sim1_Ndem,2); sim1_Nrep=mean(aux_sim1_Nrep,2);
act_sharedem=Ndem./(Ndem+Nrep);
sim_sharedem=nanmean(aux_sim_Ndem./(aux_sim_Ndem+aux_sim_Nrep),2);
sim1_sharedem=nanmean(aux_sim1_Ndem./(aux_sim1_Ndem+aux_sim1_Nrep),2);
[sim_f,sim_xi] = ksdensity(sim_sharedem);
[sim1_f,sim1_xi] = ksdensity(sim1_sharedem);
[act_f,act_xi] = ksdensity(act_sharedem);
plot(act_xi,act_f,'b'); hold on; plot(sim_xi,sim_f,'r'); plot(sim1_xi,sim1_f,'g');
